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Abstract We introduce a simple model of hard spheres with gravitational interactions, for which 
we study a suitable scaling limit. Usual extensive properties are maintained notwithstanding the long 
range of gravitational interaction. We show that a local thermalization spontaneously emerges within a 
microcanonical description of the stationary state. In the considered scaling limit, the resulting density 
profile can be determined in a hydrostatic approach. 
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1 Introduction 

It is widely accepted that, in order to even attempt a statistical mechanics treatment of a system of 
TV masses that mutually interact via Newtonian forces in three dimensions, one has to confine them 
inside a box of volume A to avoid evaporation. Confinement in a volume is a natural requirement 
also for a gas of particles interacting via short-range forces, but looks inappropriate for a model that 
should describe the Universe, a system that occupies full space. However, this seems a necessity for 
any thermodynamical theory of a self-gravitating system in the absence of cosmological expansion pQ . 

Moreover, a regularization of the interaction at short distances must be introduced to prevent 
the divergence of thermodynamic potentials [21G2II1I3I1] • In the case in which the Newtonian potential 
is softened at short distances and masses have no spatial extension, the minimal energy scales like 
TV 2 and H-stability is violated [7J. A mean-field (continuum) limit can, however, be performed by 
scaling the coupling constant as 1/-/V, which yields an extensive energy 18 . It should be remarked that 
energy remains non additive and, as a consequence, one finds that different statistical ensembles give 
inequivalent predictions {91ITU1I5] . leading to interesting phenomena like negative specific heat in the 
microcanonical ensemble |11| [2"] . If particles have instead a finite radius and cannot therefore penetrate 
each other, the minimal energy scales like N 5 ^ 3 in three dimensions, violating extensivity. 

In this paper we consider a self-gravitating system confined in a box and made of equal mass, equal 
shape and equal volume particles. Particles cannot penetrate each other as a result of a hard core 
short-range repulsion among them. As a consequence, collisions among the particles are elastic and, 
therefore, no aggregation phenomenon takes place. Particles collide elastically also with the walls of 
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the confining box. Being then the dynamics energy and particle number conserving, it is appropriate 
to consider the system in the microcanonical ensemble. The microcanonical distribution in phase space 
is a stationary solution of Liouville equation, so that it indeed describes a stationary state of the 
system B The hard core interaction prevents the divergence of the entropy. However, because of the 
long-range nature of the gravitational interaction, energy is a priori non extensive and should again 
scale like N 5 ^ 3 . We will introduce a specific scaling limit in which we recover the extensivity of the 
energy. In that scaling limit, it turns out that fluctuations of the potential energy are small compared 
to the average when N — > oo. Moreover, we show that if the energy per particle is large enough, local 
thcrmalization occurs due to hard sphere collisions. The one-body mass distribution is then found to 
obey a Boltzmann-like formula in terms of the mean-field gravitational potential, where the average 
kinetic energy plays the role of temperature. The fact that a local equilibrium is established, combined 
with a suitable separation of the length scales associated with hard core and gravitational interactions 
respectively, allows us to justify the introduction of a hydrostatic equation which express the balance 
between hard-sphere pressure and gravitational attraction. 

We stress that the present scaling provides a finite mass density, while the ratio of the kinetic 
energy to the gravitational energy can take arbitrary values. Therefore, it should be well suited for 
describing a wide class of physical situations. Notice that other types of scaling have been considered in 
the literature [TWT] , but they describe more specific situations like that of an infinitely diluted system 
for instance. Furthermore, once the scaling limit has been taken, the hard-core cutoff can be removed 
providing a meaningful limit as it would be desirable. In fact at large enough energy per particle we 
recover the well-known model of the isothermal gravitational sphere pTl[T7l[T8] . The situation is quite 
different from that studied in Ref. [3] for self-gravitating hard spheres in contact with a heat-bath in 
the canonical ensemble : when the size of the hard spheres is sent to zero while the other parameters 
are kept fixed, the equilibrium state becomes a Gaussian in velocity and a Dirac <5-peak in position, in 
both two and three dimensions. Here, in our scaling limit, collapse is avoided because the gravitational 
interaction energy between two spheres at contact vanishes. 

The paper is organized as follows. Section [5] presents the model and discuss the appropriate scaling 
continuous limit. The scaling properties of the potential energy are presented in Sec.[3j followed by the 
discussion of the mass distribution and the fluctuation of the potential energy. Section 0] presents how 
a Boltzmann-like formula for the inhomogencous density emerges spontaneously, emphasizing a local 
thcrmalization in the microcanonical ensemble. The hydrostatic approach is then discussed in Sec. [SJ 
We finally present our conclusions and draw some perspectives in Sec. [6] 



2 Definitions 

2.1 Hard spheres with gravitational interactions in the microcanonical ensemble 

We consider a classical gravitational model made of N identical hard spheres with mass m and diameter 
a, enclosed in a spherical box of volume A — 47ri? 3 /3. The corresponding Hamiltonian reads 

N 2 i 

with the two-body interaction potential 

v(r) = oo for r < a , v(r) = —Gm 2 /r for r > a . (2) 

We consider that the previous system is isolated and does not exchange energy with some thermo- 
stat. Thus, its energy is fixed and equal to some value E. We assume that the corresponding stationary 

Use of the microcanonical distribution amounts to assume, roughly speaki ng, the hypothesis of molecular 
chaos leading to some ergodicity property. In fact, as argued in Refs. |13U14) . gravitational systems may be 
trapped in metastable states through their dynamical evolution, while the corresponding transition times 
grow exponentially fast with respect to the number of particles. That mechan ism might induce a breakdown of 
ergodicity, as suggested by results obtained from numerical simulations 15 , 16 . Therefore, as far as applications 
to real astrophysical systems are concerned, the microcanonical distribution has to be handled with some care. 
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state is described within the microcanonical ensemble. The corresponding distribution of positions and 
momenta of the N particles in the canonical phase space reads 

/micro(ri,—.rjV)Pi,—»PJv) = A N 5(E-H N ) , (3) 
while the total number of microstates is 

Q[E, N,A)=A N [ TT d 3 ri d 3 Pi <J(£ - H N ) , (4) 

JA N xR3 N i 

where An is some normalization constant, which is not relevant for our purpose. 

The microcanonical distribution /micro remains stationary under time evolution generated by Hamil- 
tonian H^. We assume that there is no other conserved quantity than energy, like the total orbital 
momentum for instance [20]. We stress that Q(E,N,A) is finite for a > 0, while it diverges for point 
particles with a = for N > 3 as shown in Ref. [3]. 



2.2 The scaling continuous limit 

Since we are interested in the properties of a system with a large number of particles, it is useful to 
consider some limit where N — > oo. For ordinary systems with short-range interactions, that limit 
is nothing but the usual thermodynamical limit where both the energy per particle E/N and the 
particle density p p = N/A are kept fixed. In that limit, the physical parameters which describe one 
particle are also kept fixed. Because of both the attractive and long-range natures of gravitational 
interaction, such a limit would provide a collapsed state with non-extensive properties. In order to 
describe other physical situations of interest, other limits have been introduced in the literature, like 
the one describing an infinitely diluted system |19) . 

Here, we want to build a scaling limit when N — > oo, which describes an infinite continuous fluid 
with the usual extensive properties. For that purpose, we consider that the parameters which define the 
particles vary with N like power laws, while the gravitational constant G is not rescaled contrarily to 
some mean- field scalings introduced in the literature [5]. Therefore, we set a = d$N a and m — mgN 6 , 
while the size R of the spherical box is chosen to diverge as ^o-^ 7 with 7 > so that the system indeed 
becomes infinitely extended. The particle density p p = N/A behaves then as iV 1 ^ 37 . Since the particle 
are hard spheres, it is essential to keep the packing fraction 77 — irp p a 3 /6 bounded. By imposing that 
77 remains constant, we find a first constraint 

1 - 3 7 + 3a = . (5) 

Moreover, we impose that the mass density p = mp p remains also constant, so a second constraint 
arises, 

1 - 3 7 + S = . (6) 

Eventually, the extensivity of the gravitational energy GM 2 /R with the total mass M = Nm, provides 
the third constraint, 

2 + 28 - 7 = 1 . (7) 

The three exponents are then readily determined as a = —2/15, 5 = —2/5 and 7 = 1/5. Thus, the 
power laws defining the required scaling limit when N — > 00, denoted SL, are 

Cm 2 

R = N 1 'H , m = N-^ 5 m , a = N' 2 ' 15 ^ , E = e (8) 

with parameters £q, tuq, do and e fixed. 

Within the above SL, particles become infinitely small and light, while their inner mass density 
remains constant. Their number density p p = (3/(47t^q))A^ 2 / 5 diverges, while the mass density p = 
mp p = 3mo/(47r£g) is indeed kept fixed, as well as the packing fraction 77 = Trp p a 3 /6 — c?§/(8^q). 
That limit clearly describes an infinite continuous medium, with an extensive total energy E. The 
corresponding stationary state is controlled by two independent dimensionless parameters, namely e 
which is the energy per particle in units of Gtoq/^o, and the packing fraction 77 = dp/ (8^). Notice 
that, within the considered scaling, the mean free path I = 1/ a 2 p p ) remains proportional to the 
diameter a of the hard spheres, I — aj (6^/2 rf). Therefore, for high dilutions such that 77 1, I is 
much larger than a. 
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3 Scaling properties 

3.1 H-stability and extensivity of potential energy 

For any allowed configuration, the potential energy 

is larger than that of the collapsed configuration where the N hard spheres make a single cluster with 
size I/coii ~ iV 1//3 cr, which is of order —Gm 2 N 2 /L co \\. In the scaling limit, the collapse radius -L co ii 
diverges as N 1 ' 5 do and this provides the classical version of H-stability 

V N > -Cns^N , (10) 
«o 

where Chs is a positive real number entirely determined by the geometrical arrangement of the hard 
spheres at maximum packing. On the another hand, the potential energy should reach its maximum 
when all particles are homogeneously distributed on the spherical surface of the box. That maximum 
behaves as the gravitational energy of an empty sphere carrying the surface mass density Nm/ (AirR 2 ), 
so we expect the upper bound 

V N <-^N. (11) 

According to bounds (fTUf and (TTTj) . Vn remains of order N for any configuration. This implies 
extensive lower and upper bounds for the average potential energy (Vn), where (A) denotes the micro- 
canonical average of any observable A weighted by the distribution ([3]). Therefore, the nonextensive 
behaviours of the average potential energy encountered in other limits, like the usual thermodynamical 
limit for instance, do not occur in the present scaling limit. Accordingly, it seems quite plausible that 
the extensity of (Vjv) is ensured here, namely 

(Vn) Cm 2 , . . . 
^ 1 # = V U(£ ' r?) ' (12) 

where u(s,rj), the average potential energy per particle in units of is a well-behaved function 

of the intensive dimensionless parameters e and rj. We stress that a rigorous mathematical derivation 
of the extensive behaviour ([12"]) is beyond our scope. 



3.2 Mass distributions 

For further purposes, it is useful to introduce the n-body mass distributions. Let consider first the 
one-body mass distribution which reads, 

p(r)=m(J2S(r i -r)) . (13) 

i 

According to the scaling properties of the potential energy described above, we can reasonably expect 
that p(r) takes a well-defined shape in the SL. That statement can be made more precise as follows. 
For dimensional and spherical symmetry reasons, we can recast p(r) as 

p(r) = p(r) = pg^ } (g; e, rj) (14) 

with q = r/R and where the function (q; £, rj) is dimensionless. Notice that an explicit expression for 
9^n (?! e ' r l) can ^ e readily obtained in terms of multiple integrals upon the dimensionless positions q.; = 
Ti/R. Now, we make the assumption that gy(q;e, rj) goes to some well-defined function g^\q;s,rj) 
when N — > oo with parameters q,E,rj fixed, so 

lim p(qR) = pg {1) (g; e, rj) . (15) 
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The physical content of that scaling property, as well as its limitations, will be discussed further. 
The two-body mass distribution is 

(r, r') = m 2 (]T *(r« - r)<5(r, - r')) , (16) 

while similar definitions hold for higher order n-body mass distributions with n > 3. The scaling 
property (|15l) can be extended to such mass-distributions, namely 

lrnip^(Rq,Rq') = p 2 g^(q,q';e,r ] ) , (17) 

and so on. Of course, the limitations evocated above also apply to the present scaling behaviours. 

Notice that the extensive behaviour (|T2"j) is consistent with the above scaling assumptions for the 
mass distributions. This is easily seen by starting from the integral expression of the average potential 
energy 

(Vn) = -\( d 3 rd 3 rV 2 >(r,r')-^- . (18) 

Z J A 2 \ r - r I 

After making the variable changes r — > i?q, r' — > i?q', and using the scaling behaviour (fTTf for 
/c/^-Rq, i?q'), we find that (Vn)/N is indeed an intensive quantity in the SL given by formula (|12p 
with 

u ^^) = -^l d 3 qd 3 qV 2) (q,q';^)^^- (19) 
327r A<i,g'<i |q-q| 



3.3 Fluctuations of potential energy 



Let us now consider the fluctuations of the potential energy around its average. Similarly to for- 
mula l|18[). such fluctuations can be rewritten as 



K) [{V N )f = \ [ d 3 rd 3 ry 2 )(r,r')-^ I + / d 3 rd 3 r'd 3 rV 3 )(r, r', r") 

pW(r,r',r",r"')-p (2) (r,rV 2) (r",r"') 



m 



d 3 rd 3 r'd 3 r"dV" 

A 4 



|r — r'||r — r" 

G 2 



J r — r' r" — f 



^2P) 



where p^(r, r', r") and p^ A \r, r', r", r'") are the three- and four-body mass distributions. In the scaling 
limit, the first two terms in expression (|20|) can be estimated by replacing p^ 2 \r, r') and p^(r, r', r") 
by their scaling behaviours in terms of g^ 2 \c{, q'; e, 77) and (q, q', q"; e, rj). The first term is then 
found to be of order N°, while the second one behaves as N. 

If we estimate the third term in expression (1201) by using only the scaling behaviours of the involved 
mass distributions, we obtain not surprisingly a iV 2 -behaviour. However such an estimation should over- 
estimate the exact behaviour, because the involved correlations [p^ (r, r', r", r'") — p^ (r, r')p^ (r", r'")l 
can reasonably be expected to become rather small compared to p 4 , for spatial configurations (r, r', r", r'") 
describing a large part of the integration domain A 4 . A precise estimation cannot be performed at this 
level, so here we only assume that the correlation contribution ((54)) should grow slower than iV 2 . 

The previous simple arguments suggest that the square of fluctuations ([2U|) should become small 
compared to N 2 in the SL, i.e. 

(V 2 )-{(V N )] 2 ~o(N 2 ). (21) 

Notice that the corresponding fluctuations for an ordinary system with short-range interactions at 
thermodynamical equilibrium, do satisfy the behaviour (I21[) since they are proportional to N. 
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3.4 Scaling decomposition ansatz 

In the following, we will have to perform averages of quantities involving the potential energy Vn- 
Taking into account the extensivity of its average, as well as the behaviour (|2ip of its fluctuations, we 
propose to decompose Vn as 

V N -> (V N ) + W N , (22) 

with (Vn) = O(N) and Wn = o(N) in the SL. We assume that such decomposition ansatz holds for 
most probable configurations which mainly determine the averages of interest, so it should provide the 
exact behaviours in the cases studied further. 

Notice that both the extensivity of (Vn) and the fluctuations behaviour ([2~Tj) . although quite plau- 
sible, are not well established at this level. In fact, we will show a posteriori that the mass distributions 
and correlations inferred from the decomposition ansatz, are such that the a priori assumptions about 
(Vn) and ((V N ) — [(Vn)] 2 ) are indeed satisfied. This provides some kind of consistency check of the 
derivations, but of course it does not constitute a proof. Moreover, in the calculations of some averages, 
subtle correlations involving Vn and other observables might occur, so the above decomposition ansatz 
would no longer work, even if (Vn) and ((V N ) — [(Vjv)] 2 ) do behave as expected. 



4 Emergence of local thermalization 

In this section, we study the inhomogeneous mass density p(r) given by average (fT3"|) for energies per 
particle e > —1/2. First, we express p(v) in terms of the gravitational potential created by all the 
particles except one which is fixed at position r. Then, in the SL, we rewrite exactly the correspond- 
ing spatial average by exploiting both the extensivity and e > —1/2. Applying the decomposition 
ansatz ((22)) to that exact expression, we obtain a Boltzmann-likc formula for p(r) in the SL, with 
some temperature which emerges naturally. We conclude by a few comments, in particular about the 
corresponding formula for the one-body distribution /W(r, p) in phase space. 



4.1 Introduction of the gravitational potential 



As a starting point, we compute the average (113[) with the microcanonical distribution /micro given by 
expression Then, the standard integration over the momenta of the N particles leads to 

N 

p{r) = B(E,N,A) / Y[d 3 i i [E-V N (T,T2,...,T N )] m ^- 1 xe{E-V N (T,T 2 ,...,T N )). (23) 

In that expression, is the usual Heaviside function such that = 1 for £ > and = 
for £ < 0, while the normalization constant B(E, N, A) ensures that the spatial integral of p(r) over 
the box does provide the total mass Nm of the system. The conditions |r,; —Fj\ > cr, arising from the 
hard sphere interaction, apply to any pair of particles, in particular those including particle 1 fixed at 
ri = r. 

For e > —1/2, it turns out that (E — Vn (r, r 2 , rjy)) is positive for any spatial configuration 
thanks to the upper bound pip for the potential energy. Then, the Heaviside function can be replaced 
by 1 in the expression (|23[). That simplification is crucial for further transformations. In particular, let 
us introduce the gravitational potential ^(r|r 2 , rjv) at r created by the (N — 1) particles located at 
r 2 , ...,rjv. According to the decomposition 

V/v(r,r 2 , -,rjv) = V N -i(r 2l -,r N ) + m^(r|r 2 , ...,rjv) , (24) 
we can exactly rewrite formula (|23[) as 



N ,- , n 37V/2-l 

3/2 - m<P 



p(r) = B(E, N, A) [ dpN-i TT 0(|r 4 - r\/a - 1) X [E - V N -x] 



1 - 



E - V N -i 



, (25) 
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where d/ijv_i denotes the unnormalized microcanonical measure for the spatial configurations (r2, rjy) 
of a system made of (N — 1) particles, enclosed in the same spherical box with volume A and the same 
energy E as the genuine system with N particles, 

N 

dM;v-i = II d$r ' II 9 ^ - V M° - 1) - Vn-i(v2, r^)] 3 ^- 1 ^ 2 " 1 . (26) 

i=2 i<j 

Notice that, if s < —1/2, then Vjv-i an d & are coupled through the constraint E — Vn-i — m<P > 0, 
so the formula ([23)) a priori no longer holds. 



4.2 Exploiting the extensivity of potential energy 

Let us consider identity 



1 - 



E - V N -i 



3JV/2-1 



= exp <^ (37V/2 - 1) In 



1 - 



E - V N -i 



(27) 



Since m<P = 0(1) and E — Vjv-i = 0(N), we can expand the logarithm in powers of m<I>/ (E — Vjv-i)- 
After multiplication by factor (3N/2 — 1), we see that the linear term provides a contribution of order 
0(1), while higher powers provide vanishing contributions when N — > oo. Accordingly, we obtain for 
any spatial configuration 



1 



m<!> 



3JV/2-1 



exp 



f 3Nm<P 1 
\ 2{E-V N - 1 ) J 



E - Vjv-i 

in the SL. Inserting that asymptotic behaviour inside the r.h.s. of formula 

iv 



(28) 



we infer in the SL 
3Nm<P 



p(r) ~ B(E, N, A) jf^ d MJV -i f[ 0(|r< - r|/a - 1) x [E - V N ^ exp|- ^|^ -y j 



(29) 



We stress that asymptotic expression (|2T)]) is exact, provided that e > —1/2. The extensivity of the 
potential energy for any spatial configuration, which follows from the particular scaling defined here, 
plays a crucial role in the derivation of that behaviour. Within other scalings, which do not preserve 
that remarkable extensivity property, formula (|29p would break down. 



4.3 Boltzmann-like formula 

If we introduce the kinetic energy K^-i of the (N — 1) particles for a given spatial configuration, 
the exponential factor on the r.h.s. of the asymptotic expression (f2"9"| . can be seen as some kind of 
Boltzmann factor. However, at this level, the corresponding temperature fluctuates. Here, we show 
that such fluctuations can be neglected by applying the fluctuation ansatz described in Section [3] 

The integral with the measure d/iN-i in the r.h.s. of formula (|29p is proportional to the micro- 
canonical average of the quantity 

^ 0(\r t -v\la- l)[E Vv-d'/a exp {- -Jg™^ } (30) 

for the system with (N — 1) particles. If we introduce the corresponding microcanonical average 

/ d/x/v-iVjv-i 
{V N -x) N -x = J -^f (31) 

J A 1 *- 1 
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of Vat-i, we can rewrite 

E-V N - 1 =E-(y ff - 1 ) If - 1 -W I f- 1 , (32) 

where Wm-i denotes the deviation of Vn—i with respect to its average for a given spatial configuration. 
Now, according to the fluctuation ansatz, we assume that for the most probable configurations which 
provide the main contributions to the average of quantity (|30[) . Wm-x remains small compared to N, 
so 

[E - V N ^ 2 = [E — {V N ^) N ^f' 2 [l + o(l)] , (33) 



and 

3N 3N 



2[E-V N -i] 2[E-(V N - 1 ) N - 1 ] 



o(l) . (34) 



Moreover, the average potential energy (Vjy_ i)jv— l for the system with (N — 1) particles differs by 
a term of order O(l) from its counterpart (Vat) for the genuine system with N particles. Thus, if we 
define the microcanonical temperature T(E, N, A) through the usual relation 

3NT 

E-{V N ) = {K„) = —-, (35) 



we eventually find in the SL 

3/2 „ N 



p(r) ~ B(E, N, A) 



3NT(E, N, A) 



/ d^ A r_ 1 xTTe»(|r l -r|/cr-f)expi - 



T{E, N, A) 



(36) 



The presence of the usual Boltzmann factor exp(— m<P/T) in the asymptotic behaviour (j3"rJ)) does 
not mean that the whole system is in contact with a thermostat. We recall that the system is isolated 
with a well-defined energy E, as illustrated by the presence in behaviour (l36l) of the microcanonical 
measure dpN-i which defines the integration over spatial configurations. However, the emergence of 
temperature T(E, N, A) can be traced back to some thermalization at a local level as argued and 
exploited in the next section. That temperature is indeed an intensive positive quantity in the SL, 
namely T(E,N,A) reduces to Gmg/^o times a dimensionless function T*(e,rj) when N — > oo while e 
and r\ are the parameters which remain fixed. The interpretation in terms of a local thermodynamical 
equilibrium is confirmed by the analysis of the one-body distribution /W(r, p) in phase space. Starting 
from the corresponding microcanonical definition, and using again both the extensivity of the potential 
energy and the fluctuation ansatz, we find that /^(r, p) can be rewritten, in the SL, as p(r) times 
the Maxwell distribution proportional to exp(— p 2 /(2mT)). 

Eventually, we point out that the asymptotic formula (|36[) relies on the fluctuation ansatz, the 
validity of which might break down for some sets of parameters (e, r\) which determine the stationary 
state of the infinite system. Formula (j3SJ) might then fail for states with large collective fluctuations, 
as discussed in the next section. 



5 Hydrostatic approach 

5.1 Emergence of local equilibrium and separation of scales 

The emergence of a local thermodynamical equilibrium at temperature T, can be interpreted as result- 
ing from the collisions between the hard spheres 0- The gravitational interactions between the particles 
inside a finite volume surrounding r, can be safely neglected in the considered scaling limit since the 
mass of each particle then vanishes. In particular, notice that the gravitational interaction energy be- 
tween two spheres at contact, —Gm 2 /a, vanishes as N~ 2 / 3 in the SL, so it becomes small compared to 
the thermal kinetic energy T which remains finite. Thus, the local equilibrium is entirely determined 
by the hard-core interactions, namely by the local packing fraction rj(r) = rjp(r)/p. 

2 Before the SL is taken, namely for N finite, the gravitational interactions between close particles inside a 
small volume surrounding r should also contribute to the thermalization process. 



9 



According to the above argument, at the local level around point r, gravitation should intervene 
only through the gravitational potential 

^(r) = -^d 3 rV(r')^| (37) 

created by the whole mass distribution p(r'), which can be viewed as some external potential for the 
corresponding local system. In the expression f|36[) , this amounts to replace the fluctuating potential <P 
by 4>{r), namely to apply some kind of fluctuation ansatz to <P by assuming that its fluctuations vanish 
in the SL. 

Within the present picture, the typical length scales become completely separated in the SL. On the 
one hand, the density p(r) is expected to vary on a length scale of order R, as made explicit through 
the scaling behaviour (JTSJ). On the other hand, the local correlation length A(r) is controlled by hard 
sphere interactions, so it takes the form [2"T1|2"2"] 

A(r) = (7&gfa(r)) (38) 

where £hs( 7 ?) is some dimensionless function of rj. In the SL with q = r/R fixed, since r/(Rq) goes to 
a finite value, A(i?q) behaves as a and vanishes as iV~ 2 / 15 . Thus, the density indeed varies on a scale 
infinitely larger than the local correlation length. 



5.2 The coupled hydrostatic and gravitational equations 

The emergence of a local thermodynamical equilibrium, combined to the separation of length scales, 
lead us to write down the hydrostatic equilibrium equation 

VP(r) = -p(r)V<Xr) , (39) 

which expresses the balance between pressure and gravitational forces. Furthermore, in that equa- 
tion, the local pressure P(r) reduces to that of an homogeneous gas of pure hard spheres without 
gravitational interactions, at temperature T and number density p(r)/m, i.e. 

P(r) = IM ms(71p{r)/p) (40 ) 

m 

where pus is the dimensionless hard-sphere pressure which depends only on the local packing frac- 
tion 7j(r) = r)p(r)/p. 

The present analysis strongly suggests that the density profile in the SL is a priori exactly deter- 
mined by the coupled equations ([3T|) and (|3"9"1) . with pressure P replaced by its hard-sphere expres- 
sion P0|) . together with the total mass constraint 

/ d 3 rp(r) = Nm , (41) 
J A 

and the total energy constraint 

s nt 1 r 

— + J^vp{v)^v)=E. (42) 

Notice that all those equations and constraints can be recast in terms of the dimensionless variable q = 
r/R and of the dimensionless density profile (q; e,rj) = limsL p(Rq) / P- We recall that temperature T 
is not a given parameter, so the corresponding dimensionless temperature T*(e, if) = T/(Gwiq/^o) has 
to be determined through the resolution of the whole system of equations where e and rj are the fixed 
control parameters. 

If the gravitational long-ranged interactions can be treated as a mean-field level through the in- 
troduction of the average potential 0(r), the hard-core short-ranged interactions play a crucial role in 
the determination of the density profile. Thus, the celebrated mean-field theory for point particles [TTl 
[TSlHTj , does not provide the exact density profile in the present SL. The corresponding equations can 
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be retrieved from the present analysis by replacing in the hydrostatic equation p9[) the hard sphere 
pressure by its ideal counterpart. After an obvious integration, this provides the Boltzmann expression 

p(r) = Cpexp{-^}, (43) 

where C is some dimcnsionless normalization constant which remains to be determined. The resulting 
mean-field coupled equations for p(r) and 0(r) can be also recast in terms of the dimensionless variable 
q = r/R and of the dimensionless density profile p(Rq)/p. Now e is the sole control parameter, and it 
can be rewritten as E / (GM 2 / R) where the total mass is M = Nm. The properties inferred from that 
so-called model of the isothermal sphere have been widely studied in the literature [23] . We stress that 
the presence of the hard-core interactions should invalidate, even at a qualitative level, some mean-field 
predictions. Indeed, since the local packing fraction rj(r) = r\p(v)j p cannot exceed [2"4"1[2"5] the maximal 
packing fraction rj ma , x — 0.7405..., the exact density p(r) is bounded everywhere by prj max /r], while 
some mean-field features are directly related to a possible unbounded increase of the core density p(0). 



5.3 Consistency checks 

The derivations which ultimately provide the coupled equations Q37p and (|39p , involve various a priori 
assumptions. For consistency purposes, it is of course necessary to check such assumptions within the 
global picture which sustains the hydrostatic approach. 

A first assumption concerns the extensivity property (IT21 of (Vjv), namely the existence of a well- 
defined potential energy per particle u(e, 77) in the SL. Within the hydrostatic approach, (Vat) reduces 
to the self-gravitational energy of a sphere with mass density p(r), 

V seli = i^d 3 rp(r)0(r) = -~ J dWr'p(r)p(r')|^^- . (44) 

Since p(r) does satisfy the scaling property (fT5")l . V sc \f is indeed extensive in the SL. However, notice 
that the exact expression (TTSf of (V/v) can be recast as 

(V N ) = v scli + v COII , (45) 

with 



V con - = -\! d 3 rdVp( 2 > T >(r, v')^- (46) 
and the truncated mass distribution or mass correlations, 

p( 2 ' T )(r,r')=p (2) (r,r')-p(r)p(r'). (47) 

If the coupled equations (|37| and (|39|) do not give access to the two-body mass correlations, the 
existence of a local thermodynamical equilibrium entirely determined by hard sphere interactions 
allow us to introduce a simple description of such correlations. Because of the hard core, p( 2 ' T ) (r, r') 
reduces to — p(r)p(r') for |r — r'| < a. Moreover, we expect a decay of p( 2 ' T )(r, r') over the hard-sphere 
local correlation length A(r). As argued above, A(i?q) becomes proportional to a in the SL. Thus, the 
contribution of the vicinity of point r = i?q to 

dVp (2 < T) (r,r')^ (48) 

is of order Gp 2 a 2 = 0(N~ 4 / 15 ) in the SL. The remaining contribution to the integral (|48p of points r' 
such that |r — r'| ^> a can be roughly estimated by replacing p( 2,T )(r, r') by a constant times mp(r)/A 
which does not depend on r'. That spread homogeneous approximation is inspired by the sum rule 



/ d 3 rV (2 ' T) (r,r') = -mp(r) , 

J A 



(49) 
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which follows from particle conservation. The corresponding contribution to integral (|48|) is then of 
order GrapR 2 / A = 0(N~ 3 ^ 5 ) which becomes small compared to that of the region |r — r'| ~ a. 
Accordingly, we find that the correlation energy (|36"j) is of order Gp 2 R 3 a 2 = 0(N 1 ^ 3 ). Thus, (Vjv) is 
indeed extensive, and the potential energy per particle is entirely given by the self part V sc u , namely 

U (^) = -J^ [ d 3 qd 3 q'3 (1) (q;£^)5 (1) (q';e^)^^. (50) 

32tt 2 J q < 1>ql <i |q-q'l 

Let us consider now the fluctuations ({V$) — [(Vjv)] 2 ) of the potential energy, which are given by 
formula (|20p . We can rewrite the first two terms as mean- field like terms 

\ _( 2 d 3 rdVp(rMr')^^ (51) 

and 

J a3 d 3 rd 3 r-d 3 r-p(r)p(r>(r-) |r _^ r m _ r//| , (52) 

plus the corresponding correlation terms associated with [p( 2 )(r,r') — p(r)p(r')] and [p^ 3 '(r,r',r") — 
p(r)p(r')p(r")}. Since p(r) does satisfy the scaling property (|T5|) . the mean-field contributions (|5T|) 
and (|52|) are of order TV and N respectively. The corresponding two- and three-body correlation con- 
tributions are readily estimated within a simple modelization of the truncated distributions analogous 
to the one used above for analyzing contribution (|4l)f to the average (Vn) itself. They are found to 
become small compared to their mean-field counterparts in the SL. 

It remains to estimate the contribution of the third term in expression (|2"0")) of the fluctuations 
({V£)-[{V N }} 2 ). If we define 

n™(r",r'") = p^(r,r',r",r"')/p^(r,r') - p {2 \r",r"') , (53) 
we can rewrite that four-body correlation term as 

\ J^r'p^M^ / A2 d 3 r"d 3 r'"nS,(r",r- )F ^ . (54) 

Similarly to the case of p ( - 2 ' T - ) (r, r'), we expect that, for two given points r and r', rr 2 ^ (r",r"') takes 
non- vanishing values for spatial configurations such that one or more relative distances |r" — r|, |r" — r'|, 
|r"' — r'|, |r"' — r"| is of order a. The corresponding contributions to the integral 

/ A / 3 r''d 3 r'%^(r'',r''' )F ^ (55) 

are then readily estimated along similar lines as above when analyzing the two-body correlation 
term (|48p. For each of the four regions when one of the points r" or r'" is close to either r or r', 
we find a contribution to integral (|55p of order Gp 2 a 3 R 2 = 0(1). After integration over r and r', the 
corresponding contribution to the term (154")) is O(N). All the other spatial configurations (r",r"') for 

which hard-sphere correlations determine rr ', (r , r ), ultimately provide contributions o(N) to ([54^1 . 
Eventually, it remains to determine the contributions of regions such that all distances |r" — r|, |r" — r'|, 
|r"' — r'|, |r"' — r"|, are large compared to a. As for the case of p^ 2,T - 1 (r, r'), we again use a spread 
homogeneous approximation inspired by the sum rule 

f d 3 r"d 3 r"'n^,(r", r'") = 2(3 - 2N)m 2 , (56) 

J A 

namely we replace n^, (r", r'") by a constant proportional to Nm 2 /A 2 . This provides a contribution 
of order 1 to integral (p)B"|) . and ultimately a contribution O(N) to (IM)) . Accordingly, the four-body 
correlation term (f54"|) is of order N, so fluctuations {(Vjy) — [(Vat)] 2 ) indeed grow slower than iV 2 , in 
agreement with the assumption (|2"Tj) . 
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5.4 Further limitations 

We stress that there are other implicit assumptions in the hydrostatic approach, which might break 
down, at least for some sets of values for the control parameters (e, r/). In particular, if the local packing 
fraction exceeds some critical value ?7lCj the system is expected to undergo a phase transition from 
a liquid to a cristalline phase [2~l~|f2"2"]. Then, the density p(r) would oscillate with a spatial period of 
order a in the bulk, for r = i?q with |q| < 1. Therefore, the hydrostatic equation which assumes 
local homogeneity on a large scale compared to a should be no longer valid [j. This might occur for 
sufficiently low values of e and not too small values of r\. On the contrary, for not too low values of e, 
and rj sufficiently small, we can expect that the local packing fraction ?7(r) never exceeds the critical 
value f?LCj so the system remains in a fluid phase everywhere as implicitly assumed in the hydrostatic 
approach. 

Another important assumption in the derivations relies on the existence of a single most probable 
state which provides the leading contributions to the averages of the quantities of interest, while the 
corresponding corrections are associated with fluctuations which remain relatively small. Within the 
inferred hydrostatic approach, this means that the solution of the coupled equations ((37)) and f|39[) 
has to exist and to be unique. Such existence and uniqueness are in fact not guaranteed. This is 
well illustrated by the analysis of the point-like version of those equations with r\ = : there are 
no solutions for e < e m i n ~ —0.335 [18] . while two and more solutions can coexist in some range 
£min < £ < £ m ax [23] . For j) 7^ 0, numerical studies [261127] using approximate simple forms of the hard- 
sphere equation of state pns (v) indicate that some of the previous features are still observed, so they are 
not specific to the point-like case 77 = 0. In fact, they are related to the attractive nature of gravitational 
interactions which favors the collapse of the system. The absence of solutions to the coupled equations 
(l37l) and (|39| might be related to the flatness of the microcanonical distribution. When two or more 
solutions coexist, one would expect some kind of multimodal shape for that distribution. In any of such 
cases, there appears macroscopic fluctuations which cannot be neglected with respect to the averages 
themselves. Thus, the whole scheme leading to the hydrostatic approach should then break down. 
A proper many-body microcanonical analysis of the system, which might then undergo some phase 
transitions [22] > remains a challenging and rather difficult problem. 

The tendency to collapse induces either the violation of the condition r)(r) < r/LC near the core 
r = 0, or the lack of the uniqueness property of the hydrostatic solution. That mechanism is of 
course related to Jeans instability [28 , which is controlled by the Jeans length ij cx {T / {Gpm)) 1 / 2 . 
Within the SL, Lj reduces to the size R of the system times [T*(e, 7?)/3] 1//2 . When e — > 00 at a fixed 
V < ?7LC) since the potential energy is always negative, the dimensionless temperature T*(e, 77) also 
diverges. More precisely, the lower and upper bounds (|10[) and (1111) on that potential energy, imply the 
behaviour T*(e,r]) ~ 2e/3 when e — > 00 at a fixed r\ < 771,0 Accordingly, the ratio Lj/R also diverges, 
so previous limitations arising from the tendency to collapse should not occur. Thus, for e sufficiently 
large and rj sufficiently small, the hydrostatic approach can be reasonably expected to provide the exact 
density profile. In fact, the gravitational interactions can then be treated as a small perturbation, so 
the density p(r) should be almost uniform and close to p, in relation with the quite plausible limit 
behaviour <jrW (q; s, rj) — !• 3/(47r) when e — > 00 with q and rj < t/lc fixed. When e is decreased, 
temperature T*(e, rj) should decrease so the ratio Lj/R should also decrease. When Lj becomes of the 
order of the size R, or even smaller than R, Jeans-like instabilities should play an important role and 
invalidate the assumptions underlying the hydrostatic approach. Therefore, there should exist some 
threshold value e c (rf) depending on 77, below which the hydrostatic approach fails. Notice that, if the 
derivation of the hydrostatic approach has been carried out with the condition e > —1/2, this does not 
mean that e c (rj) is necessarily smaller than —1/2, as suggested by the point-like case 77 = for which 
no mean-field solutions exist for e < e m i n ~ —0.335. 



3 Notice that close to the wall of the spherical box, i.e. for r almost equal to R, the mass density p(r) varies 
on a scale a. The corresponding shape does not show in the SL of p(qR) for fixed q < 1. 
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6 Conclusion 

We have presented in this paper a novel approach to the statistical mechanics of self-gravitating parti- 
cles in the microcanonical ensemble which emphasizes that gravitational interactions can be treated at 
the mean-field level. The main idea is that, while the gravitational constant is not rescaled, an appro- 
priate scaling of the parameters of the particle allows us to define a non-trivial equilibrum state with 
the usual H-stability and extensivity properties. If this scaling is introduced through simple physical 
requirements, it turns out that, rather unexpectedly, it also ensures the spontaneous emergence of the 
local thermodynamical equilibrium (|43[) . 

Several comments and open questions are however important to draw at this stage. 

i) We show here that, locally, the gas of self-gravitating particles is thermalized through collisions, 
so the present scaling defines a temperature. Moreover this finding justifies a posteriori some previous 
attempts to study self-gravitating particles within a canonical ensemble [21[S]. However, that thermal- 
ization is expected to occur only at sufficiently large enough energy per particle. Furthermore, although 
the thermalization is recovered at the local level, there is still no clear definition of an external ther- 
mostat at the macroscopic level, with which the system is exchanging heat or kinetic energy. Indeed, 
as it is well known, the energy being conserved, a gain (or loss) of kinetic energy can be obtained by a 
decrease (or increase) of potential energy. It is still debated nowadays [351 how one can use a canonical 
framework with a clear physical meaning. 

ii) The hard-core regularization is essential, on the one hand for avoiding the collapse of the finite 
system [5J, and on the other hand for ensuring through collisions the emergence of a local equilibrium 
in the scaling limit. Once that limit has been taken, hard-core effects intervene in the local pressure 
which ultimately determines the density profile, at least for a large enough energy per particle. At this 
level, the hard core can be removed and we recover the density profile for point particles determined 
within the familiar mean-field description of the isothermal sphere. 

iii) We have considered throughout the paper a large enough energy per particle e, but of course 
states with negative e are of strong interest. Indeed, we have discovered that positive values prevent 
the Jeans length to be small enough and thus hinder any equilibrium collapsed state. A non vanishing 
fraction of matter will stay on the boundary of the system, even if the number of particles N diverges 
to infinity. Statistically, a large part of these particles will have momenta pointing outwards, leading 
to the partial evaporation of the system. This effect is not described within the present model in which 
all particles are maintained artificially into the spherical domain. To avoid this problem, one might 
consider negative energies, which would lead to sufficiently small Jeans length to constrain particles to 
stay far from the boundaries: the fraction of evaporating particles might therefore be vanishingly small. 
Interestingly, this effect is reminiscent of two particles in a Keplerian system, which are restricted to 
confined elliptic trajectories only for negative energies. 

Both long-range and short-range effects, as well as confinement, can therefore be traced back in the 
results derived within this new scaling approach to statistical mechanics in the microcanonical ensem- 
ble. They are at the origin of the difficulties, but also of the interests of this fascinating problem. We 
are however convinced that this new scaling approach is an important step forward for understanding 
the complete statistical mechanics treatment encompassing the possibility to tackle phase transitions 
and fragmentations of a system of self-gravitating particles. 
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